%==============================================================================
% This code computes equilibrium with and without CBDC under different CBDC 
% as the economy gets cashless, i.e., the department stores that accept
% both cash and electronic payment moves online and accept only electronic
% payment
% Cournot competition in both deposit and loan markets
% with a perfectly competitive interbank market and
% banks free entry.
%
%  COPY RIGHT: Yu Zhu
%==============================================================================



clear;close all
options = optimset('maxiter',5e5,'tolX',1e-10,'tolfun',1e-10,'maxfunevals',5e5);

%==========================================
%Define parameters and functions
%==========================================
name = 'LN87_08_final'
load(['pars_post_' name '.mat'])
name = [name,'_entry'];
com_opt = 1 % loan market competition 1 if cournot 0 if perfect competition

location='figures/'

alpha1=pars(4);alpha2=pars(5);alpha3=pars(6);beta=pars(7);epsilon=0.001;sigma=pars(1);
cost=pars(8);                   
B=pars(2);
A=pars(10);  %productivity
eta=pars(9);                                       %bank's cost to manage deposits
theta=pars(11);
Nb=pars(3); %number of banks
rq_ratio=pars(12);
i_reserve = pars(13);
year_option = pars(14); % 1: using banking data from 2007 and 2008. 0: using banking data from 2014-2019
%% Define the name of the saved files: last two number indicates reserve requirement ratio and the interest on reserves
if year_option==1
    pi_m = 0.03346;pi_h=pi_m;
    name = [name,'_', num2str(floor(rq_ratio*100)),'_',num2str(floor(i_reserve*100)) '_07_08'];
else
    pi_m = 0.01515;pi_h=pi_m;
    name = [name,'_', num2str(floor(rq_ratio*100)),'_',num2str(floor(i_reserve*100))];
end

picbdc=pi_m;

pi_eff = (1+pi_m)/(1+i_reserve)-1;

nrhogrid =10000;
i_m = 1/beta*(1+pi_m)-1;i_h = 1/beta*(1+pi_h)-1;

%================================================================
%Define functions
%================================================================





u = @(x) ((x+epsilon).^(1-sigma)-epsilon.^(1-sigma))./(1-sigma); %DM utility
du = @(x)(x+epsilon).^(-sigma); %derivative of DM utility
d2u=@(x)(x+epsilon).^(-sigma-1)*(-sigma);
c=@(x)x;                                                         % DM cost
dc=@(x)1;%   derivative of DM cost
d2c=@(x)0;
y_star = fminsearch(@(x)(du(x)-dc(x)).^2,1);                     % y_star
g_func=@(x) (1-theta).*u(x)+theta.*c(x);                         %trading proctocol 
dg_func=@(x) (1-theta).*du(x)+theta.*dc(x);
d2g_func=@(x)(1-theta).*d2u(x)+theta.*d2c(x);
y_grid = 0:0.001:y_star;
p_grid = g_func(y_grid);
p_star = p_grid(end);
y_func = @(p)interp1(p_grid,y_grid,min(p,max(p_grid)));          %trading quantity 
p_func=@(x)min(x,max(p_grid));                                   %trading payment
lambda_func = @(z) max(du(y_func(z))./dg_func(y_func(z))-1,0);          % liquidity premium

f=@(k)A*k.^(eta);%production
if com_opt==1
    df_inv = @(rho,Nb)(A*eta*(1-(1-eta)/Nb)./(1+rho)).^(1/(1-eta));
    df=@(k,Nb)A*eta*(1-(1-eta)/Nb)*k.^(eta-1); %marginal production
else
    df_inv = @(rho,Nb)(A*eta./(1+rho)).^(1/(1-eta));
    df=@(k,Nb)A*eta*k.^(eta-1); %marginal production
end


dlambda_func=@(z)((d2u(y_func(z)).*dg_func(y_func(z))-du(y_func(z)).*d2g_func(y_func(z)))./dg_func(y_func(z)).^3).*(z<=p_grid(end));
dinv=@(z)beta*(alpha2*lambda_func(z)+1);
ddinv=@(z)beta*alpha2*dlambda_func(z);
z_grid = linspace(0,p_grid(end),100000);
lambda_grid = lambda_func(z_grid);
gamma_grid=linspace(0,20,nrhogrid);
x0=[0,0,0];
x0_R=x0;
opt=1;
icbdc_grid=linspace(-0.01,0.05,1000);
D_grid = linspace(0,max(p_grid),5e3);
pass_variables
ygrid = D_inv_demand(i_m,D_grid,func_var);
alpha2_grid = alpha2 + linspace(0,0.03,100);
alpha3_grid = alpha3 - linspace(0,0.03,100);
eq_picbdc=(1+picbdc)/(1+0)-1;


%====================================================================
%Compute entry cost
%====================================================================
 y=SS_eq_noCBDC(func_var,pi_eff,Nb,ygrid);

psi= beta*(alpha2*lambda_func(y(2))+alpha3*lambda_func(y(1)+y(2)))+beta;
drate=1/psi-1;
loan= fzero(@(x)df(x,Nb)-1-y(end),[1e-10,100]);

eq_final_nocbdc1=[y,drate,loan];

 y=SS_eq_noCBDC(func_var,pi_eff,Nb+1,ygrid);

psi= beta*(alpha2*lambda_func(y(2))+alpha3*lambda_func(y(1)+y(2)))+beta;
drate=1/psi-1;
loan= fzero(@(x)df(x,Nb+1)-1-y(end),[1e-10,100]);

eq_final_nocbdc2=[y,drate,loan];
rl_rate1 = (1+eq_final_nocbdc1(:,3))/(1-(1-eta)/Nb)-1;
rl_rate2 = (1+eq_final_nocbdc2(:,3))/(1-(1-eta)/(Nb+1))-1;
welfare_banker1 = max(1+rl_rate1-1/(1+pi_eff),0).*eq_final_nocbdc1(:,end)-eq_final_nocbdc1(:,2)...
    -(cost-(1/(1+pi_eff))).*eq_final_nocbdc1(:,2)./(1+eq_final_nocbdc1(:,4));

welfare_banker2 = max(1+rl_rate2-1/(1+pi_eff),0).*eq_final_nocbdc2(:,end)-eq_final_nocbdc2(:,2)...
    -(cost-(1/(1+pi_eff))).*eq_final_nocbdc2(:,2)./(1+eq_final_nocbdc2(:,4));

welfare_banker1=welfare_banker1/Nb;
welfare_banker2=welfare_banker2/(Nb+1);
entry_cost = (welfare_banker1+welfare_banker2)/2; %entry cost defined as the middle point of the profits with Nb and Nb+1 banks.


for i=1:length(alpha2_grid)
    i
    alpha2=alpha2_grid(i);
    alpha3 = alpha3_grid(i);
    % Compute equilibrium without CBDC
    pass_variables
    for nb = Nb:1:200
    if i==1
        y=SS_eq_noCBDC(func_var,pi_eff,nb,ygrid);
    else
        y=SS_eq_noCBDC_fast(func_var,pi_eff,nb,x0,opt);
    end
    
    i_m=(1+pi_m)/beta-1;
    psi= beta*(alpha2*lambda_func(y(2))+alpha3*lambda_func(y(1)+y(2)))+beta;
    drate(i)=1/psi-1;
    loan= fzero(@(x)df(x,nb)-1-y(end),[1e-10,100]);
    
    eq_final_nocbdc(i,:)=[y,drate(i),loan];
    rl_rate1 = (1+eq_final_nocbdc(i,3))/(1-(1-eta)/nb)-1;
    profit = max(1+rl_rate1-1/(1+pi_eff),0).*eq_final_nocbdc(i,end)-eq_final_nocbdc(i,2)...
    -(cost-(1/(1+pi_eff))).*eq_final_nocbdc(i,2)./(1+eq_final_nocbdc(i,4));
    if profit/nb<entry_cost
        eq_final_nocbdc(i,:) = eq_cand;
        profit_nocbdc(i,:) = [profit_cand,nb-1];
        psi=psi_cand;
        break;
    end
    eq_cand = eq_final_nocbdc(i,:);
    profit_cand = profit;
    psi_cand = psi;

    end
    x0=y;
   
    %===============================================
    %Compute equilibrium with CBDC
    %===============================================
    icbdc_real  = (1+picbdc)/(0+1)/beta-1;
    if icbdc_real>psi/beta-1
        eq_cbdc(i,:)=eq_final_nocbdc(i,:);
        drate_cbdc(i)=eq_final_nocbdc(i,4);
        eliquidity(i,1) = eq_final_nocbdc(i,2);
        profit_cbdc(i,:) = profit_nocbdc(i,:);
    else
        tic
        z_eq= SS_eq( alpha1,alpha2,alpha3,z_grid,lambda_grid,i_m,icbdc_real);
        
        lsupply_max=(1-rq_ratio)*(1+icbdc_real)*beta*z_eq(1,2);
        rho_hat=max(1/(1+pi_eff)-1, ((0+1)/(1+picbdc)-rq_ratio/(1+pi_eff)+cost)/(1-rq_ratio)-1);
        for nbb = profit_nocbdc(i,2):-1:1
            ld_low= fzero(@(x)df(x,nbb)-1-rho_hat,[1e-10,100]);
            drate_cbdc(i)=(0+1)/(1+picbdc)-1;
            toc
            if ld_low<lsupply_max
                deposit =  ld_low/(icbdc_real+1)/beta/(1-rq_ratio);
                eq_cbdc(i,:)=[z_eq(1),deposit,rho_hat,drate_cbdc(i),ld_low];
            else
                rhocbdc = df(lsupply_max,nbb)-1;
                
                eq_cbdc(i,:)=[z_eq,rhocbdc,drate_cbdc(i),lsupply_max];
                
            end
            eliquidity(i,1) = z_eq(2);

            rl_rate1 = (1+eq_cbdc(i,3))/(1-(1-eta)/nbb)-1;
             profit_cbdc_nb = max(1+rl_rate1-1/(1+pi_eff),0).*eq_cbdc(i,end)-eq_cbdc(i,2)...
                -(cost-(1/(1+pi_eff))).*eq_cbdc(i,2)./(1+eq_cbdc(i,4));
            if profit_cbdc_nb/nbb >entry_cost
                profit_cbdc(i,:) = [profit_cbdc_nb,nbb];
                break;
            end
        end

    end
    
    
    
    
end
Output_noCBDC = alpha1*eq_final_nocbdc(:,1)+alpha2_grid(:).*eq_final_nocbdc(:,2)+alpha3_grid(:).*min(eq_final_nocbdc(:,1)+eq_final_nocbdc(:,2),p_star)...
    +B+eq_final_nocbdc(:,end)+f(eq_final_nocbdc(:,end))-eq_final_nocbdc(:,2)...
    +(eq_final_nocbdc(:,2)./(1+eq_final_nocbdc(:,4))-eq_final_nocbdc(:,end))*(1+i_reserve)/(1+pi_m);
Output_CBDC = alpha1*eq_cbdc(:,1)+alpha2_grid(:).*eliquidity+alpha3_grid(:).*min(eq_cbdc(:,1)+eliquidity,p_star)+B+eq_cbdc(:,end)+f(eq_cbdc(:,end))-eq_cbdc(:,2)...
    +(eq_cbdc(:,2)./(1+eq_cbdc(:,4))-eq_cbdc(:,end))*(1+i_reserve)/(1+pi_m);
output_index=Output_noCBDC/Output_noCBDC(1);
output_indexCBDC=Output_CBDC/Output_CBDC(1);




%% plot graph for cashless limit 


macro_implication=figure
psi_D = eq_final_nocbdc(:,2)./(1+eq_final_nocbdc(:,4));
psi_D_R = eq_cbdc(:,2)./(1+eq_cbdc(:,4));
Delta_grid = alpha2_grid-alpha2_grid(1);
subplot(3,3,1)
hold on
h1=plot(Delta_grid/alpha3_grid(1)*100,(psi_D./psi_D(1)-1)*100,'linewidth',2);
h2=plot(Delta_grid/alpha3_grid(1)*100,(psi_D_R./psi_D_R(1)-1)*100,'--red','linewidth',2);

xlim([0,5.5]) 
xlabel('\Delta')
ylabel('Deposit')

subplot(3,3,2)
hold on
plot(Delta_grid/alpha3_grid(1)*100,(eq_final_nocbdc(:,end)./eq_final_nocbdc(1,end)-1)*100,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,(eq_cbdc(:,end)./eq_final_nocbdc(1,end)-1)*100,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Loan')

subplot(3,3,3)
hold on
plot(Delta_grid/alpha3_grid(1)*100,(output_index-1)*100,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,(output_indexCBDC-1)*100,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Output')

ndrate=((1+eq_final_nocbdc(:,4))*(1+pi_m)-1)*100;
ndrateR=((1+eq_cbdc(:,4))*(1+pi_m)-1)*100;
subplot(3,3,4)
hold on
plot(Delta_grid/alpha3_grid(1)*100,((1+eq_final_nocbdc(:,4))*(1+pi_m)-1)*100,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,((1+eq_cbdc(:,4))*(1+pi_m)-1)*100,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Deposit Rate')

if com_opt==1
    rl_rate = (1+eq_final_nocbdc(:,3))./(1-(1-eta)./profit_nocbdc(:,2))-1;
    rl_rateR = (1+eq_cbdc(:,3))./(1-(1-eta)./profit_cbdc(:,2))-1;
else
    rl_rate = eq_final_nocbdc(:,3);
    rl_rateR = eq_cbdc(:,3);
end

nlrate=((1+rl_rate)*(1+pi_m)-1)*100;
nlrateR=((1+rl_rateR)*(1+pi_m)-1)*100;

nlrate_interbank = ((1+eq_final_nocbdc(:,3))*(1+pi_m)-1)*100;
nlrate_interbankR = ((1+eq_cbdc(:,3))*(1+pi_m)-1)*100;


subplot(3,3,5)
hold on
plot(Delta_grid/alpha3_grid(1)*100,nlrate,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,nlrateR,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Loan Rate')

subplot(3,3,6)
hold on
plot(Delta_grid/alpha3_grid(1)*100,nlrate-ndrate,'linewidth',2)
plot(Delta_grid/alpha3_grid(1)*100,nlrateR-ndrateR,'--red','linewidth',2)
xlim([0,5.5]) 
xlabel('\Delta') 
ylabel('Spread')

hL = subplot(3,3,8);
poshL = get(hL,'position');     % Getting its position

lgd = legend(hL,[h1,h2],'No CBDC','With CBDC', 'NumColumns',2);
set(lgd,'position',poshL);      % Adjusting legend's position
axis(hL,'off');                 % Turning its axis off

saveas(macro_implication,[location,'cashless_', name ],'epsc')
saveas(macro_implication,[location,'cashless_', name],'png')

id = find(eq_cbdc(:,end)./eq_final_nocbdc(:,end)>1,1,'first');

fraction = Delta_grid(id)/alpha3_grid(1);

disp(['A zero-interest CBDC promotes bank intermediation if ' num2str(fraction*100) '% Type 3 transactions reject cash'])
